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We study the time evolution of occupation numbers for interacting Fermi-particles in the situation 
when exact compound states are "chaotic". This situation is generic for highly excited many- 
particles states in heavy nuclei, complex atoms, quantum dots, spin systems and quantum computer 
models. Numerical data show perfect agreement with a simple theory for the onset of thermalization 
in close systems of interacting particles. 
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It is known that highly excited states can be treated 
as "chaotic" ones in many-body systems, such as com- 
plex atoms |jj , multicharged ions , nuclei Q, spin sys- 
tems |^|,|5| and quantum computer models |^[7j. This 
happens due to a very high density of many-particle 
states which strongly increases with an increase of en- 
ergy. For example, in the case of n Fermi-particles occu- 
pying the finite number m of "orbitals" (single-particle 
states), the total number N of many-body states grows 
exponentially fast with an increase of number of parti- 
cles, N = m\jn\[m — n)\ ~ exp(co«). Correspondingly, 
the density pf of those many-body states which are di- 
rectly coupled by a two-body interaction, also grows very 
fast. Therefore, even a relatively weak interaction be- 
tween the particles can lead to a strong mixing between 
unperturbed many-body states ("basis states"). As a re- 
sult, an exact (perturbed) eigenstate is represented by a 
chaotic superposition of a large number of components of 
basis states 



10 . The number of principal basis com- 



ponents in such chaotic eigenstates can be estimated as 
Np C ~ T/D where T is the spreading width of a typi- 
cal component that can be expressed through the Fermi 
golden rule, and D is the mean interval between the 
many-body levels. 

In this paper we extend the quantum chaos approach 
to the problem of a time evolution of an initially excited 
basis state. This initial state may contain one or sev- 
eral excited electrons above other electrons in the ground 
state of a quantum dot or atom. 

Exact many-body eigenstates \k) of a Hamiltonian 
H = Ho + V of interacting Fermi-particles can be ex- 
pressed in terms of simple shell-model basis states | /) of 
H , 



\k) 



E c fl/> 
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h 



7> 
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(1) 



Here |0) is the ground state, af is the creation operator 
and (7^ are components of an exact chaotic eigenstate 
k) in the unperturbed basis, that is formed by a residual 
interaction V. 

Below we consider the time evolution of the system, 
assuming that initially (t = 0) the system is in a specific 
basis state \i). This state can be expressed as a sum over 
exact eigenstates, 



E^ (fc) i fc >> 



(2) 



therefore, the time-dependent wave function reads as 

*(*) = E C i k) °f k ) I/) ex P H£ (fe) i). (3) 
kj 

Here £ (fc) are the eigenvalues corresponding to the eigen- 
states \k >. The sum is taken over the eigenstates \k > 
and basis states |/ > (hereafter we assume % = 1). 

The occupation number n a of a single-particles state 
a for the wave function (||) of the system is given by the 
expression, 

n a = (*(t)\n a Mt)) = E (<l\na\q) (s<*> + S<*>) . (4 ) 



with n a = aJ^Qa. 



Here Sq d ^ is the diagonal term, 



5 f = Elcf ) | 2 |cf| 2 ) 



(5) 



and Sq stands for off-diagonal terms, 



<?(/') = ]T c[ h) C^c[ p) C^ exp [i - £W) t 



k^ P 



(6) 
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Since at t — the only unperturbed basis state \i) is 

excited, the occupation numbers n a {t) are equal to or 

1 (for Fermi systems). For large times, t — > oo, the terms 

( fi) 

in the off-diagonal sum S q oscillate very rapidly as a 
function of spacings EW-E<p), therefore, S q fl) tends to 
zero. Thus, the asymptotic distribution n a (oo) of occu- 
pation numbers is determined by the diagonal term only. 

If the number of principal components in the wave 
function is large, one can replace Sq® by its average 
value, 

Sf> = F{Ei,E q ) tt^F(Ei,EW)F{E q ,EW) (7) 

k 

w J F(Ei,E)F(E q ,E)p(E)dE. (8) 

Here p{E) = D^ 1 (E) is the density of energy spectrum, 
and 



F(E h E^) = \C, 



(fc)i 



(9) 



is the so-called F -function that characterizes the shape 
of eigenstates in the unperturbed basis. This F-function 
has been studied in details for different models (see, e. 
g. flP^Jl^,p| and the review @). As was shown in 
P,|l5j, in many cases the F-function can be approximated 
by the following expression (see also fL6f), 



F(E l ,E)p(E)=B 



exp 



(E-Ec) 




(10) 



In this expression B is the normalizing constant that can 
be found from the relation J F(Ei, E)p(E)dE = 1, and 
Ti is the spreading width of a compound state \i). The 
energy E c stands for the center of effective energy band 
a of the Hamiltonian matrix, and may not coincide with 

If residual interaction V is not very strong, the spread- 
ing width is determined by the standard golden rule, 
Ti rs 2n\Vif\ 2 pf -C a. Here pf is the density of final 
basis states |/) that are coupled to \i) by the interaction 
V. In the limit T <C o the parameter a 2 is equal to the 
variance a 2 of pf. In the opposite limit r>ir the value 
of a 2 is approximately the same as the variance A 2 ^ of 
the strength function F(E h E)p(E), see Ref. @. Note 
also, that a 2 is always finite due to a finite range of an 
interaction in the energy representation. 

The form of F{Ei,E) in this limit r <C a has simple 
Breit-Wigner (BW) shape (apart from long tails that are 
highly non- universal Q|), 



F{E h E) = 



2tt P (E - E t ) 2 +r 2 /4' 



(11) 



where r = Tj. In this case the function F(Ei,E q ) that 
gives the final shape of the time-dependent wave function, 



has also the BW form Qrj with the width r = 2I\ Note 
that fluctuations of the diagonal term are small, 
s (d) _ s (d) _ S ( p/y/N^ c . Thus, for a large number of 
principal components N pc the occupation numbers both 
in compound eigenstates and in large-time asymptotic of 
the wave function, have the self-average property. As a 
result, we obtain the following expression for the asymp- 
totic values of occupation numbers, 



i a (oo) = 22(q\h a \q)F(Ei,E q ). 



(12) 



In fact, this distribution is very close to the distribution 
of occupation numbers in compound eigenstates of a sys- 
tem, 



^2{q\n a \q)F(E q ,E). 



(13) 



The properties of the distribution (|T^) have been studied 
in details HgHj^HQ for different models. 

In many cases (the two-body random interaction 
(TBRI) model fTM2|, s-d nuclear shell model p|, multi- 



ply charged ions M) the ^-distribution (13) is very close 
to the standard Fermi-Dirac distribution for the occu- 
pation numbers with a certain temperature T(E). In 
these cases the asymptotic distribution ( |l2| ) of the oc- 
cupation numbers in the time-dependent problem is also 
given by the Fermi-Dirac distribution. In other cases, like 
Ce atom p"7| where the residual interaction between par- 
ticles on some orbitals is larger than the distance between 
single-particle energy levels, the F-distribution (|l^) de- 
viates very strongly from the Fermi-Dirac distribution. 
This happens since the mean-field approximation is not 
a "good" one. However, even in this case the distribu- 
tion of n a is close to that obtained from the canonical 
distribution that takes into account a repulsion between 
electrons (see details in |l7j ) . This canonical distribution 
is also characterized by some temperature T{E). Thus, 
it is natural to term the time evolution of n a (t) as the 
process of thermalization in a closed system. 

The natural question is: how fast is this termalization? 
To answer this question, one needs to study the evolution 
of a many-body system taking into account a two-body 
nature of the interaction between particles. For this, it 
is convenient to consider the so-called cascade model (see 
details in p9[ ), that is based on the representation of 
unperturbed many-body states in the form of "classes" . 
In this picture, the first class contains those N\ basis 
states that are directly coupled by the two-body inter- 
action Vi f . Correspondingly, the second class consists of 
N2 basis states that are coupled with the initial state 
in the second order of the interaction V (therefore, the 
coupling is proportional to VikVkf), etc. 

The analytical expression for time dependence W n (t) 
of a population of the class n has been obtained in Ref. 
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p9| . The population in the "zero" class is just a proba- 
bility Wo(t) of the system to remain in an initial state. 
In the BW case, F <C a, this probability for large times 
is Wo(t) ~ exp(— Ft). In other limiting case, r > a, 
when the form of the F'-function is close to the Gaus- 
sian, the time dependence Wo(t) is more complicated 
p6[ . For example, for r > u, it was found to be 
Wo(t) = exp(— A|.i 2 ), up to a very long time. In gen- 
eral case it changes from Wo(t) = exp(— A E t 2 ) for small 
times, to Wo(t) — Cexp(— Ft) for large times and sat- 
urates near Wo(t) w 3/N pc (here C > 1 is the time- 
independent coefficient, see details in fll6[). 

Let us now consider the thermalization process. Ini- 
tially, the occupation numbers n a (0) are equal to or 1. 
In the BW regime, the first class is populated during the 
time r ~ 1/T |L9| (or, r ~ 1/Ae in the Gaussian regime, 
for r > <j). At that time the occupation numbers n a al- 
ready strongly deviate from their initial values since the 
two-body interaction can move any of two particles to 
new single-particles levels characterized by the energies 
e a . Thus, the characteristic time for a thermalization of 
the occupation numbers is determined by the population 
time r for the first class. Note, that this time is also 
a characteristic time for the "decay" of the probability 
Wo(t). To compare with, the population of all n c classes 
in a particular system requires a longer time t Uc ~ n c r 
fL9f . Thus, in the case of n c ^> 1 (e.g. in a mesoscopic 
system) the thermalization of the occupation numbers 
may occur on a smaller time scale than the onset of a 
complete statistical equilibrium. 

This suggests a simple derivation of the time depen- 
dence n a (t) for occupation numbers. From the normal- 
ization condition ^s=o W s = 1 one can find the popu- 
lation of all classes with s ^ as 2s=i = 1 — Wo- 
Now we assume that the thermalization of the occupa- 
tion numbers occurs on the time scale r. This assumption 
leads to a simple expression 

n a (t) = n a {0)W (t) + n Q (co) (1 - W (t)) . (14) 

Here n a (0) are initial occupation numbers (0 or 1), and 
n a (oo) determines equilibrium occupation numbers (e.g., 
the Fermi-Dirac distribution). 

In Figs. 1-2 we compare numerical data for n a (t) with 
the above estimate (|l4j). Two situations are studied, the 
first (Fig.l) corresponds to the BW regime (weak inter- 
action regime, r <C a) , and in the second (Fig. 2) the 
form of the F- function is close to the Gaussian (r > a). 
For numerical simulation, we have used the model with 
random two-body interaction (see details, in Jl2|) with 
n = 6 Fermi-particles that occupy m = 12 single-particle 
levels. Since the analytical expression for Wo(t) in gen- 
eral case is quite complicated, in numerical simulations 
we use the exact value of Wo(t) computed numerically. 

Overall, there is an excellent agreement between the 
theory and numerical data (apart from fluctuations that 



are neglected in the theory) . In order to simplify the pic- 
ture, we have taken the initial basis state \i) just in the 
middle of the many-particle spectrum that consists of 924 
levels. In this case all final values of the occupation num- 
bers are equal, n a (oo) = n/m =1/2 that corresponds to 
the infinite temperature T. 




FIG. 1. Time dependence of occupation numbers n a (t) for 
the BW regime (weak interaction ). Thin curves with dots 
present numerical data, thick smooth curves correspond to 
the analytical expression (|l^), see the text. Computations are 
made for the model with random two-body interaction, with 
n = 6,m = 12,7? =< V 2 /D% >« 0.003, F w 0.50, A E « 1.16 
(here Do is the mean spacing between single-particle ener- 
gies). 




FIG. 2. Time dependence of occupation numbers n a (t) for 
the case when the strength function has the gaussian form: 
77 =< V 2 /Dl >« 0.083, T w 10.5, A E « 5.8. 

We would like to point out on an interesting difference 
between two cases of a weak and strong interaction be- 
tween the particles. Specifically, for a weak interaction 
(Fig.l) the transition to equilibrium values of n a has a 
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character of damped oscillations. The number of princi- 
pal components N pc ~ Tp" 1 in this case is not very large, 
this is why there are considerable fluctuations in n a (t) 
even in the equilibrium. To compare with, the case of a 
strong interaction (Fig. 2) shows fast and monotonic tran- 
sition to thermal values of n a (oo) with relatively small 
fluctuations, see also p9[ . 

Our results demonstrate that actually there are two 
time scales in the onset of thermalization. The first one 
is determined by r. It characterizes an "initial thermal- 
ization", and allows one to use Eq.(|l4|) for a descrip- 
tion of the time dependence of occupation numbers. For 
larger times, damped quantum oscillations (with period 
T ~ n c r) may occur in the transition to a complete equi- 
librium [GM9I. 



In conclusion, we present a simple theory for the onset 
of thermalization in closed systems of interacting parti- 
cles. The theory allows to describe the time dependence 
of the occupation numbers, and shows the relaxation to 
an equilibrium distribution. Numerical data obtained for 
the model of a random two-body interaction between fi- 
nite number of Fermi-particles, demonstrate a very good 
correspondence to the theoretical predictions. Our re- 
sults can be used in different applications, such as com- 
plex atoms, heavy nuclei, atomic clusters, quantum dots, 
etc. 
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